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We obtain a new value for the QCD coupling constant by combining lattice QCD simulations 
with experimental data for hadron masses. Our lattice analysis is the first to: 1) include vacuum 
polarization effects from all three light-quark flavors (using MILC configurations); 2) include third- 
order terms in perturbation theory; 3) systematically estimate fourth and higher-order terms; 4) 
use an unambiguous lattice spacing; and 5) use an (!J(a^)-accurate QCD action. We use 28 different 
(but related) short-distance quantities to obtain a^^{Mz) ~ 0.1170(12). 

PACS numbers: 11.15.Ha,12.38.Aw,12.38.Gc 


An accurate value for the coupling constant as in quan¬ 
tum chromodynamics (QCD) is important both for high- 
energy phenomenology, and as an input for possible theo¬ 
ries beyond the Standard Model. Numerical simulations 
of QCD using lattice techniques, when combined with ex¬ 
perimental data for hadron masses, have provided some 
of the most accurate values for the coupling constant P| . 
The precision of these determinations has been limited, 
however, by two factors. One was our inability to include 
the effects of realistic light-quark vacuum polarization in 
QCD simulations. The other limitation was the lack of 
third and higher order terms in the perturbative expan¬ 
sions used to extract as- In this paper we present the 
first lattice QCD determination of the coupling constant 
that includes realistic vacuum polarization effects from 
all three light quarks, and perturbation theory through 
third order, with systematic estimates of fourth order 
and beyond. Consequently, our final results are, by far, 
the most accurate from lattice QCD and among the most 
accurate from any method. This work uses gluon config¬ 
urations from the MILC collaboration , and builds on 
a joint effort by several groups 0. 

Effects from light-quark vacuum polarization are quan¬ 
titatively important, but also very costly to simulate. 
Previous simulations included contributions from only 
u and d quarks, no s quarks, and used quark masses that 
were 10 times too large or larger. Our analysis includes 
effects from all three light quarks, with much smaller 
u and d masses — so small that our results become ef¬ 
fectively mass-independent. This is possible because of 
a new discretization of the light-quark action . Heavy- 
quark polarization is negligible and is ignored here 0 . 

The Lorentz-noninvariant ultraviolet regulator greatly 
complicates high-order perturbation theory in lattice 
QCD. To manage this complexity, we automated the 
generation of Feynman integrands using computers, and 


evaluated the Feynman integrals numerically on large- 
scale parallel computers. These techniques allowed us to 
evaluate perturbative coefficients through third order . 

To extract the coupling constant from our lattice QCD 
simulation, we (with our collaborators) first tuned the 
theory’s five parameters to reproduce experiment for five 
well-measured quantities; the details are in We 

used lattices that were approximately 2.5 fm on a side 
with three different lattice spacings a, where a~^ was 
1.144(31), 1.596(30), and 2.258(32) GeV. The s quark 
masses we used for each of the three lattice spacings were 
0.082/a, 0.05/a, and 0.031/a, respectively. We used u 
and d masses as small as ms/5, except on the coarsest 
lattice where we used rris/lO. The gluon configurations 
were produced using an improved gluon action and the 
new light-quark action. Our coupling-constant analysis 
is the first to use 0(a^)-accurate actions. 

The lattice spacing a is one of the five simulation pa¬ 
rameters, and the most important in our analysis because 
it sets the simulation’s mass scale. In our earlier as anal¬ 
yses, we set the lattice spacing by comparing a simulated 
T mass splitting {e.g., T' —T) with experiment. Here we 
continue this practice, but, for the first time, the lattice 
spacings derived from our T splitting have been shown 
to agree with spacings derived from a wide variety of 
other physical quantities: ten in all, including the pion 
and kaon leptonic decay constants, the Bs mass, and the 
Q baryon mass 0, . All of these different quantities 

agree on the lattice spacing to within 1.5-3%. 

Having an accurately tuned simulation of QCD, we 
used it to compute nonperturbative values for a variety 
of short-distance quantities, each of which has a pertur¬ 
bative expansion of the form 
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Y = '^Cnay{d/a) (1) 
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where c„ and d are dimensionless a-independent con¬ 
stants, and av{d/a) is the (running) QCD coupling con¬ 
stant, for 71/ = 3 flavors, in the V scheme |3,l3- Given 
the coefficients c„, we determine av(d/a) such that the 
perturbative formula for Y reproduces the nonperturba- 
tive value from the simulation. 

We computed Cn for n < 3 using Feynman diagrams. 
We estimated higher-order coefficients by simultaneously 
fitting results from different lattice spacings to the same 
perturbative formula. This is possible because the cou¬ 
pling avid/a) changes value with different lattice spac¬ 
ings a: 

= -Pooty - Piay - - Pa^v -( 2 ) 

where the Pi are constants Q . We parameterize the run¬ 
ning coupling in our fits by its value at a specific scale — 
ao = av(7.5GeV)—and integrate Eq. numerically 
to obtain values at other scales. 

We used a constrained fitting procedure, based upon 
Bayesian methods, for our fitsj^. Given simulation re¬ 
sults Yi ± avi for three different lattice spacings ai ± Uai, 
we minimized an augmented function. 
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by varying ao, and the c„ and at. Terms after the first 
in are “priors” that constrain the fit parameters to a 
reasonable range. The fits explored values of c„, for ex¬ 
ample, centered around cT with a range specified by crc„ ■ 
For n < 3, we set clT to the values obtained from our nu¬ 
merical evaluations of the relevant Feynman diagrams, 
with (Tc„ equal to the uncertainties in those evaluations. 
We set = 0 for 4 < n < 10, and CTc„ = CTc where CTc 
was determined using the Empirical Bayes procedure de¬ 
scribed in 0. (Typically this procedure set somewhat 
larger than the optimal value found for |c 4 |.) We ignored 
terms with n > 10. The prior constraint on the coupling 
constant, ao, was 0.20to;io, or equivalently —1.6±0.7 for 
log(ao); it had negligible impact on the fits. 

The simplest short-distance quantities to simulate are 
vacuum expectation values of Wilson loop operators: 

= i (0| ReTrPe-*®^"-^'^" |0), (4) 

where P denotes path ordering, Afj, is the QCD vector po¬ 
tential, and the integral is over a closed ma x na rectangu¬ 
lar path. Wilson loops are perturbative when ma and na 
are small. We computed perturbative coefficients for six 
small loops , “measured” them nonperturbatively in 
simulations with each of our three lattice spacings, and 


did fits to perturbation theory for each loop. We also 
evaluated Wilson loops for two non-planar paths 0 : 



The fits revealed that high-order coefficients in the 
Wilson loop expansions are larger than we expected: for 
example, we find 


logWii = -3.068 ay(3.33/a) (1 - 1.068 ay (6) 
-f 1.69(4) a^ - 5(2) a^ - 1(6) a^. ■■■); 
logWi 2 = -5.551 ay (3.00/a) (1 - 0.858 ay (7) 
-f 1.72(4) a^ - 5(2) a^ - 1(6) a^ ■ ■ ■) . 

The large 5ay corrections are needed if perturbation the¬ 
ory is to agree with simulation results for all three lattice 
spacings. The coupling ay(3.33/a) ranges between 0.21 
and 0.29 for our lattice spacings, so 5ay is 5-12% of the 
full result. Each Wilson loop we examined had correc¬ 
tions of order this size. 

These large coefficients reduce the accuracy of our final 
results. There are two ways to reduce the size of these 
coefficients. One is to “tadpole improve” Wmn by divid¬ 
ing by where uq = (ITii)^^^ 0- The other is to 

examine Greutz ratios of the loops rather than the loops 
themselves 0 . Each procedure significantly reduces the 
high-order coefficients we obtain when we refit to results 
from our three lattice spacings: for example. 


log|^-^^ = 0.949 ay(1.82/a) (1-^0.160(2) ay (8) 

-0.54(8) a^-2(l)a^-0(2) a^ •••); 

log(^) = -1.323 ay(1.21/o) (1-0.39(1) ay (9) 
V ^22 J 

-b0.3(2) air - 2(1) air + 0(2) a^, ■■■). 


These expansions are typical of the 7 tadpole-improved 
loops and 6 Greutz ratios that we examined. Each 
has smaller ay coefficients, but also significantly smaller 
scales for the ays. Over our range of lattice spacings, 
ay(1.21/a), for example, ranges between 0.33 and 0.68, 
and therefore 2ay is 7-60% of the final result depending 
upon the lattice spacing. Consequently results from these 
quantities are not significantly more accurate than those 
from Wilson loops. Results from the coarsest lattices, 
with large ays, carry the least weight in our fits. 

We also examined the static-quark potential, which is 
perturbative at short distances. The continuum potential 
has a particularly simple form in the V scheme: 


V(r) 


-Cf 


a\%(0.5614/r) 

r 





( 10 ) 


where Cp = 4/3 and /So = 11 — 2nfl2>. On the lattice 
we examined the quantity Vir) — V(a) since lattice arti¬ 
facts cancel almost completely in the difference. We com¬ 
puted and removed the small residual lattice artifacts in 
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FIG. 1: Values for the 5-flavor at the Z mass from each 
short-distance quantity. The dashed lines indicate our final 
result, 0.1170(12) (x^ per data point is 0.77). 


V{r) — V (a) through second order in perturbation the¬ 
ory, and fit the resulting potential with the continuum 
formula; higher-order lattice artifacts are negligible here. 
Continuum perturbation theory for V(r) — V (a) becomes 
nonanalytic, however, in fourth order, with the appear¬ 
ance of terms proportional to ay log(av) 113 . log(av) 
is small for our range of ays, so we see no evidence of 
it in our fits. Nevertheless, the presence of such terms 
suggests that our results from the potential may not be 
as reliable as those from our other quantities. We limited 
our analysis to r < 3a as otherwise the ay scales become 
too small (< 1/a). For the same reason, we discarded 
results for the potential from the coarsest lattice. 



d/a (GeV) 


FIG. 2: Values for av versus d/a (Eq. 0) from each short- 
distance quantity at each lattice spacing, with (top) and with¬ 
out (bottom) light-quark vacuum polarization. The dashed 
lines show predictions from Eq. 0 assuming av (7.5 GeV) is 
0.2082(40) and 0.1645(14) for n/ = 3 and 0, respectively. 


Finally, we extracted the coupling directly from the 
tadpole-improved bare lattice coupling, aiat/bFii, which, 
like Wilson loops, has large fourth-order coefficients. 

We extracted values for the 3-flavor coupling, ao = 
Q;v(7.5GeV), from fits to each of our 28 short-distance 
quantities. To facilitate comparison with other determi¬ 
nations, we converted our results from the V scheme to 
the MS scheme [3 , added c and b quark vacuum polariza¬ 
tion (oerturbative lv [ll| . using quark masses of 1.25(10) 
and 4.25(15) GeV |l2|). and evolved to the Z mass. The 
results from the different quantities are shown in Fig.^ 
While they are derived from the Wilson loops, 
our Greutz ratios and tadpole-improved loops provide 
coupling-constant information that is largely indepen¬ 
dent of that coming from the loops. This is because the 
highly ultraviolet contributions that dominate the loops 
largely cancel in the other quantities, making the lat¬ 
ter far more infrared (c./., {d/a)s for loops and ratios). 
Our 28 separate determinations of the scale parameter 
probe a wide range of different length scales, have very 
different sensitivities to potential nonperturbative errors, 
and, as we have discussed, have very different perturba¬ 
tive expansions. The agreement, to within our errors, 
of all 28 determinations is strong evidence that we have 
correctly identified and controlled the various systematic 
errors that could have affected our analysis. 

The weighted average of our 28 determinations gives a 
composite result of 

= 0-1170(12), (11) 

or, equivalently, 

(7.5 GeV) = 0.2082(40). (12) 

Our error estimate here is that of a typical entry in the 
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log Wii 

log W 13 /W 22 

V(V2a) - V{a) 

a ^ 

0.0007 

0.0010 

0.0010 

Cl ... Ca 

0.0001 

0.0004 

0.0004 

c„ for n > 4 

0.0008 

0.0005 

0.0004 

U ^ MS ^ Mz 

0.0001 

0.0001 

0.0001 

condensate 

0.0002 

0.0001 

0.0001 

ruu, md, ms 

0.0004 

0.0001 

0.0001 

rric, mb 

0.0002 

0.0002 

0.0002 

simulation errors 

0.0000 

0.0000 

0.0002 

total uncertainty 

0.0012 

0.0012 

0.0012 


TABLE I: Sources of the uncertainties in our final determina¬ 
tions of the coupling q^(Mz) = 0.1170(12). 

plot; combining our results does not reduce errors be¬ 
cause most of the uncertainty in each result is systematic. 

Our composite value for the coupling constant agrees 
well with the current Particle Data Group world average 
value of 0.1187(20) , but is somewhat more accurate. 

It also agrees within errors with our previous results 0 , 
and with the preliminary result of the present analysis 
presented in [Sj, where we quoted larger uncertainties 
because we lacked values for the C 3 S. 

Realistic vacuum polarization is critical to our result. 
Redoing our simulations and analysis, but with no light- 
quark vacuum polarization, gives a coupling of 0.0900(4) 
rather than 0.1170(12). Evolving to Mz increases the 
difference between 0 and 3 light-quark flavors, but the 
couplings are still lOu apart at 7.5 GeV. This is evident in 
Fig .[ 21 where we plot the av{d/a)s extracted from our fits, 
with and without vacuum polarization, for each quantity 
and for each of our lattice spacings separately, but using 
the c„s from our simultaneous fits to all lattice spacings. 

The various sources of uncertainty for different quanti¬ 
ties are elaborated in Tabled The dominant errors come 
from three sources. First is the uncertainty in the inverse 
lattice spacing a“^, which includes both statistical and 
estimated finite-a errors in the simulated upsilon splitting 
(0.5-2% depending upon a ^3)- Second are residual un¬ 
certainties in the parameters ci... C 3 from the numerical 
calculation of these coefficients. The final large source 
of uncertainty is due to uncertainties in the coefficients 
beyond third order. This error is greatly reduced be¬ 
cause we fit simultaneously to three lattice spacings; fit¬ 
ting with just a single lattice spacing, as is usually done, 
would give errors 2-5 times larger. We allowed for pos¬ 
sible effects from nonperturbative gluon and quark con¬ 
densates 0 , but these are negligible. Uncertainties in 
the c and b masses, and Monte Garlo simulation errors in 
the loop values are negligible. We corrected for the errors 
in M, d and s sea-quark masses by redoing our entire anal¬ 
ysis (loops and lattice spacings) with larger masses and 


extrapolating linearly. This is also negligible; we include 
an uncertainty in Tabled equal to the correction. 

Our coupling-constant analysis uses the most realistic 
QGD simulation to date, with, for the first time, vacuum 
polarization contributions from all three light-quarks, 
quark and gluon actions corrected through 0{a^), and 
extensive evidence that both the heavy- and light-quark 
sectors of the theory have been accurately simulated . 
It is the first to use not only third-order accurate per¬ 
turbation theory, but also systematic estimates of fourth 
order and higher. Our final results come from 28 differ¬ 
ent short-distance quantities, covering almost an order of 
magnitude in energy scales. The agreement between our 
results and the current world average demonstrates that 
the QGD of confinement is the same theory as the QGD 
of jets; lattice QGD is full QGD, encompassing both its 
perturbative and nonperturbative aspects. 
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